Calculate functional diversity metrics
Casey O’Hara
3/5/22
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
fig.height = 4, fig.width = 7)
library(terra)
library(sf)
library(data.table)
library(tidyverse)
library(here)
library(oharac)
source(here('common_fxns.R'))Summary
Using functional entities defined in a prior script based on a suite of categorical traits, in combination with spatial distribution of species from AquaMaps and IUCN, calculate functional diversity metrics based on Mouillot et al (2014) and code from Sebastien Villeger: functional richness, functional vulnerability, functional redundancy, and functional over-redundancy.
Methods
Summarize functional entity membership per cell across IUCN and AquaMaps
Included species are:
- coded for functional entity membership
- mapped in one (or both) of the species distribution mapsets
- included the trait-based vulnerability set
Inspect bin membership for FE traits
spp_info <- assemble_spp_info_df() %>%
select(species, taxon, fe_id, mob, trp, len, wcol) %>%
distinct() %>%
mutate(mob = factor(mob, levels = c('ses', 'sed', 'mob', 'mig')),
trp = factor(trp, levels = c('pri', 'sec', 'ter', 'apex')),
len = factor(len, levels = c('tiny', 'vsm', 'sm', 'med', 'ml', 'lg', 'vlg', 'huge')))
n_fes_possible <- n_distinct(spp_info$mob) *
n_distinct(spp_info$trp) *
n_distinct(spp_info$len) *
n_distinct(spp_info$wcol)
n_fes_realized <- spp_info$fe_id %>% n_distinct()Overall, with this set of traits used to define functional entities, we could have up to 512 FEs, only 339 (66.2%) of which are realized in this set of species. Note that the Mouillot 2014 paper includes two additional traits bringing the total number of possible FEs up to 5670, only 646 of which are realized by the 6316 tropical fish spp in their study.
Adult body length
knitr::kable(t(table(spp_info$len)))| tiny | vsm | sm | med | ml | lg | vlg | huge |
|---|---|---|---|---|---|---|---|
| 2174 | 5959 | 4179 | 3595 | 2555 | 1370 | 906 | 406 |
Trophic level
knitr::kable(t(table(spp_info$trp)))| pri | sec | ter | apex |
|---|---|---|---|
| 1647 | 4247 | 12810 | 2440 |
Adult mobility
knitr::kable(t(table(spp_info$mob)))| ses | sed | mob | mig |
|---|---|---|---|
| 2539 | 4378 | 10709 | 3518 |
Water col position
knitr::kable(t(table(spp_info$wcol)))| ben | bp | pel | rf |
|---|---|---|---|
| 12702 | 1130 | 2233 | 5079 |
Load species ranges and summarize to functional entities per cell
Here we use the IUCN species range maps rasterized to 10 km
Mollweide, and AquaMaps species distributions reprojected to 10 km
Mollweide. Priority: IUCN species maps, then AquaMaps with occur_cells ≥
10. These map filenames are already set up in the
assemble_spp_info_df() function.
Read in functional entity assignments and join to spatial ranges. Per cell per FE, calculate the number of species for later calculation into functional richness, redundancy, vulnerability. Break the process into chunks to avoid overloading memory.
process_spp_to_fe <- function(fe, spp_maps) { ### fe <- fe_vec[1]
dt_out <- spp_maps %>%
filter(fe_id == fe) %>%
data.table() %>%
.[ , .(n_spp_fe = length(unique(species))),
by = .(cell_id, fe_id)]
return(dt_out)
}fe_summary_file <- here_anx('func_entities', 'fe_cell_summary.csv')
# unlink(fe_summary_file)
if(!file.exists(fe_summary_file)) {
spp_info_df <- assemble_spp_info_df(fe_only = TRUE, vuln_only = TRUE)
### table(spp_info_df %>% select(species, len) %>% distinct() %>% .$len)
# huge lg med ml sm tiny vlg vsm
# 406 1370 3595 2555 4179 2174 906 5959
### Chunk out and save smaller files to avoid loading in EVERY spp map
### all at once!
n_chunks <- 40
spp_to_chunks <- spp_info_df %>%
select(species) %>%
distinct() %>%
mutate(chunk = rep(1:n_chunks, length.out = n()))
spp_clean <- spp_info_df %>%
select(species, fe_id, map_f, taxon) %>%
distinct() %>%
left_join(spp_to_chunks, by = 'species')
spp_fes <- spp_clean %>%
select(species, fe_id) %>%
distinct()
tmp_chunk_filestem <- here('tmp/fe_summary_chunk_%s.csv')
for(i in 1:n_chunks) { ### i <- 1
tmp_chunk_f <- sprintf(tmp_chunk_filestem, i)
if(file.exists(tmp_chunk_f)) {
message('Temp file ', basename(tmp_chunk_f), ' exists... skipping!')
next()
}
message('Assembling species maps and binding FEs for chunk ', i, ' of ', n_chunks, '...')
spp_chunk <- spp_clean %>% filter(chunk == i)
chunk_spp_maps <- collect_spp_rangemaps(spp_vec = spp_chunk$species,
file_vec = spp_chunk$map_f,
idcol = 'species',
) %>%
dt_join(spp_fes, by = 'species', type = 'left') %>%
distinct()
fe_vec <- chunk_spp_maps$fe_id %>% unique() %>% sort()
message('... summarizing FR map for chunk ', i, ': ', length(fe_vec), ' FEs in ',
nrow(chunk_spp_maps), ' cell observations...')
chunk_fe_sumlist <- parallel::mclapply(fe_vec, mc.cores = 20,
FUN = process_spp_to_fe,
spp_maps = chunk_spp_maps)
if(check_tryerror(chunk_fe_sumlist)) {
stop('Try-errors detected in FE summary loop!')
}
message('... binding chunk and writing chunk summary to ', basename(tmp_chunk_f), '...')
chunk_fe_summary <- chunk_fe_sumlist %>% data.table::rbindlist()
fwrite(chunk_fe_summary, tmp_chunk_f)
}
message('... reading temp chunks, binding, and summarizing by fe and cell...')
fe_tmp_fs <- list.files(dirname(tmp_chunk_filestem),
pattern = 'fe_summary_chunk', full.names = TRUE)
fe_summary_list <- parallel::mclapply(fe_tmp_fs, FUN = fread, mc.cores = 10)
fe_summary <- rbindlist(fe_summary_list) %>%
.[ , .(n_spp_fe = sum(n_spp_fe)),
by = .(cell_id, fe_id)]
message('... writing summary to ', basename(fe_summary_file), '...')
fwrite(fe_summary, fe_summary_file)
### unlink(fe_tmp_fs)
}Calculate functional diversity metrics
Calculate functional diversity metrics from IUCN and AquaMaps summaries of functional entity membership per cell.
This set of calculations replicates the results of Villeger’s function but without resorting to a presence/absence matrix, instead just relying on the original AquaMaps species-cell table.
- spp_fe_df joined to am_spp_cell;
- group by cell and functional group;
- calculate number of spp per functional group;
- group by cell;
- calculate functional entity metrics -
- Species richness
- Functional entity richness
- Functional vulnerability (% of FEs represented by a single spp)
- Weighted functional vulnerability (new: vuln of an FE is \(\frac{1}{2^{n-1}}\), take average over all FE)
- Functional redundancy (mean spp count per FE)
- Functional over-redundancy (sum of spp per FE above mean, divided by total number of spp)
fe_metrics_file <- here_anx('func_entities/fe_metrics_per_cell.csv')
### unlink(fe_metrics_file)
if(!file.exists(fe_metrics_file)) {
# system.time({
fe_sum_total <- data.table::fread(fe_summary_file)
### there are 6.5 million cells. Chop into 500k instances and mclapply to summarize
# cell_id_vec <- 1:ncell(raster::raster(here('_spatial/ocean_area_mol.tif')))
chunk_size <- 250000
n_chunks <- 6.5e6 / chunk_size
nspp_tmp_f <- here('tmp/nspp_df_tmp.csv')
if(!file.exists(nspp_tmp_f)) {
message('Calculating species richness map...')
nspp_list <- parallel::mclapply(1:n_chunks, mc.cores = ceiling(n_chunks / 4),
FUN = function(n) { ### n <- 1
message('... spp richness chunk ', n, ' of ', n_chunks, '...')
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
nspp_sum <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
data.table() %>%
.[ , .(n_spp = sum(n_spp_fe, na.rm = TRUE),
n_fe = n_distinct(fe_id)),
by = 'cell_id']
})
if(check_tryerror(nspp_list)) stop('Try errors detected!')
nspp_df <- nspp_list %>%
data.table::rbindlist()
write_csv(nspp_df, nspp_tmp_f)
} else {
message('File exists: ', basename(nspp_tmp_f))
nspp_df <- data.table::fread(nspp_tmp_f)
}
fvuln_tmp_f <- here('tmp/fvuln_df_tmp.csv')
if(!file.exists(fvuln_tmp_f)) {
message('Calculating functional vulnerability map...')
fv_list <- parallel::mclapply(1:n_chunks, mc.cores = ceiling(n_chunks / 4),
FUN = function(n) { ### n <- 1
message('... funct vuln chunk ', n, ' of ', n_chunks, '...')
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
f_vuln <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
oharac::dt_join(nspp_df, by = 'cell_id', type = 'left') %>%
data.table() %>%
.[ , .(f_vuln = sum(n_spp_fe == 1) / n_fe,
f_wvuln = mean(0.5^(n_spp_fe - 1)),
f_red = mean(n_spp_fe)),
by = .(cell_id, n_fe, n_spp)]
})
if(check_tryerror(fv_list)) stop('Try errors detected!')
f_vuln_df <- fv_list %>%
data.table::rbindlist()
write_csv(f_vuln_df, fvuln_tmp_f)
} else {
message('File exists: ', basename(fvuln_tmp_f))
f_vuln_df <- data.table::fread(fvuln_tmp_f)
}
fred_tmp_f <- here('tmp/f_red_df_tmp.csv')
if(!file.exists(fred_tmp_f)) {
message('Calculating functional overredundancy map...')
f_overred_list <- parallel::mclapply(1:n_chunks, mc.cores = ceiling(n_chunks / 4),
FUN = function(n) { ### n <- 1
message('... funct overredundancy chunk ', n, ' of ', n_chunks, '...')
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
f_overred <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
oharac::dt_join(f_vuln_df, by = 'cell_id', type = 'left') %>%
data.table() %>%
.[ , f_over := ifelse(n_spp_fe > f_red, n_spp_fe - f_red, 0)] %>%
.[ , .(f_overred = sum(f_over) / sum(n_spp_fe)),
by = 'cell_id']
})
if(check_tryerror(f_overred_list)) stop('Try errors detected!')
f_overred_df <- f_overred_list %>%
data.table::rbindlist()
write_csv(f_overred_df, fred_tmp_f)
} else {
message('File exists: ', basename(fred_tmp_f))
f_overred_df <- data.table::fread(fred_tmp_f)
}
message('Combining functional diversity metrics maps...')
fe_metrics_df <- f_vuln_df %>%
oharac::dt_join(f_overred_df, by = 'cell_id', type = 'full')
message('Combined set: ', nrow(fe_metrics_df), ' rows!')
write_csv(fe_metrics_df, fe_metrics_file)
}Map functional diversity metrics
Create and save rasters of functional diversity metrics
fe_metrics_df <- data.table::fread(fe_metrics_file)
### note that while there are 3,684,273 cells in this dataframe,
### 4273 of those are the Caspian Sea, so when those are dropped,
### it results in 3,680,000 exactly... ugh, weird
n_spp_rast <- map_to_mol(fe_metrics_df, which = 'n_spp')
n_fe_rast <- map_to_mol(fe_metrics_df, which = 'n_fe')
f_vuln_rast <- map_to_mol(fe_metrics_df, which = 'f_vuln')
f_wvuln_rast <- map_to_mol(fe_metrics_df, which = 'f_wvuln')
f_red_rast <- map_to_mol(fe_metrics_df, which = 'f_red')
f_overred_rast <- map_to_mol(fe_metrics_df, which = 'f_overred')
writeRaster(n_spp_rast, here('_output/func_entities/n_spp.tif'),
overwrite = TRUE)
writeRaster(n_fe_rast, here('_output/func_entities/n_fe.tif'),
overwrite = TRUE)
writeRaster(f_vuln_rast, here('_output/func_entities/f_vuln.tif'),
overwrite = TRUE)
writeRaster(f_wvuln_rast, here('_output/func_entities/f_wvuln.tif'),
overwrite = TRUE)
writeRaster(f_red_rast, here('_output/func_entities/f_red.tif'),
overwrite = TRUE)
writeRaster(f_overred_rast, here('_output/func_entities/f_overred.tif'),
overwrite = TRUE)Examine distributions of each raster.
Plot each raster. NOTE: These no longer include spp that are missing
from vulnerability calculations - so the maps in
_output/nspp_maps (vuln \(\cap\) FE) should be identical to these in
output/func_entities.
